# connect the GRN to the (expressed & connected) interface TF objects from Metacore

source("../scripts/hairballFunctions.r")

# arguments
args=commandArgs(TRUE)
cutoff = as.double(args[1])

# read the contextualized GRN (1 solution)
GRN = read.table(sprintf("interactions_least_pruned-cutoff%s.txt",cutoff),sep="\t",stringsAsFactors=FALSE)

# read the interactions among all possible interface TFs and TFs
int_to_TFs = read.table("../data/MetacoreIntTFobject-GRN_TFs-interactions.txt",sep="\t",stringsAsFactors=FALSE,check.names=FALSE,quote="")

# check which interface TFs are connected to pathway source nodes through expressed path
components = read.table("../data/TransReg_objectTOgeneSymbol.txt",sep="\t",stringsAsFactors=FALSE,header=TRUE)
connected_table = read.table(sprintf("table-reachableTFs-cutoff%s.txt",cutoff),sep="\t")
TFs = rownames(connected_table)
components = components[components[,1] %in% TFs,]
library(igraph)
graph_str = readRDS("../data/curated_graph.Rdata")
components = components[components[,1] %in% V(graph_str)$name,]

# take the interface objects that are expressed
pexpr = read.table("bool_MetacoreObjects.txt",stringsAsFactors=FALSE,row.names=1,sep="\t",quote="")
pe = pexpr[unique(components[,1]),2]
expr_nodes = unique(components[,1])[pe>cutoff]

# check if the GRN TFs are among the interface nodes: for such interface nodes I require all TF components DE (AND expressed)
inGRN = aggregate(components[,2] %in% c(GRN[,1],GRN[,3]), list(components[,1]), all)
DE_interf = intersect(inGRN[which(inGRN$x),'Group.1'],expr_nodes) #they need to be expressed and DE

# add interactions from non DE,expressed nodes, to the TFs already in the GRN
nodes_to_add = setdiff(expr_nodes,DE_interf)
int_to_TFs = int_to_TFs[(int_to_TFs[,1] %in% nodes_to_add) & (int_to_TFs[,3] %in% c(GRN[,1],GRN[,3])),] 
newGRN = rbind(GRN,int_to_TFs)
newGRN = unique(newGRN)

# takes interactions (only +/-) and gives adj_matrix
adj_matrix = simpleTintToAdjmat(newGRN) # see hairballFunctions.r

# prepare the pheno data
phenored = read.table(sprintf("phenored-DEbool-cutoff%s.txt",cutoff),row.names=1) # read phenotypes of contextualized GRN
to_sort = cbind(colnames(adj_matrix),phenored[match(colnames(adj_matrix),rownames(phenored)),1],phenored[match(colnames(adj_matrix),rownames(phenored)),2])
to_sort[which(is.na(to_sort))] = 2
sorted = to_sort[order(to_sort[,2]),] # phenotypes for the GRN+interface TFs network; 2=to perturb in both directions
write.table(sorted,sprintf("pheno-contextualized-cutoff%s-interface.txt",as.character(cutoff)),sep="\t",quote=FALSE,row.names=FALSE,col.names=FALSE)

# prepare the adjacency matrix
adj_matrix = adj_matrix[match(sorted[,1],rownames(adj_matrix)),match(sorted[,1],colnames(adj_matrix))]
write.table(adj_matrix,sprintf("adjmat-contextualized-cutoff%s-interface.txt",as.character(cutoff)),sep="\t",quote=F)

# list of interface nodes to perturb
nodes=unique(int_to_TFs[,1])
single_TFs = unique(components[components[,1] %in% DE_interf,2] )
to_perturb = c(single_TFs,nodes)
write.table(cbind(to_perturb),sprintf("to-perturb-cutoff%s.txt",as.character(cutoff)),sep="\t",row.names=FALSE,col.names=FALSE,quote=FALSE)

# dictionary for interface nodes that were already in the GRN
symb = components[components[,1] %in% DE_interf,2:1]
write.table(symb,sprintf("singleTF_to_node-cutoff%s.txt",as.character(cutoff)),row.names=FALSE,col.names=FALSE,sep="\t",quote=FALSE)

# define interface TFs
interfaceTFs = to_perturb
if (nrow(symb)>0){
	interfaceTFs = c(interfaceTFs[!interfaceTFs %in% symb[,1]], symb[,2])	
}
interfaceTFs = unique(interfaceTFs)
write.table(cbind(interfaceTFs),sprintf("interfaceTFs-cutoff%s.txt",cutoff),row.names=FALSE,col.names=FALSE,quote=FALSE)

# define perturbation of 3 or 4 TFs based on the number of interface TFs
n_to_perturb = length(to_perturb)
print(sprintf("Number of interface TFs to perturb: %i",n_to_perturb))
if (n_to_perturb >100) {
	print("Combinations of 1-3 TFs will be perturbed")
} else {
	print("Combinations of 1-4 TFs will be perturbed")
}

## Prepare pert.m file
pert = sprintf("pert-cutoff%s.m",cutoff)
unlink(pert)
system(paste0("touch ",pert))
l1 = paste0("\"subfix='cutoff",cutoff,"'",";\"")
system(paste0("echo ",l1," >>",pert))
l2 = paste0("\"addpath(", "'../scripts/'" ,");\"")
system(paste0("echo ",l2," >>",pert))
system(paste0("echo perturb_interface >>",pert))
system(paste0("echo perturb_interface_2nodes >>",pert))
system(paste0("echo perturb_interface_3nodes >>",pert))
if (n_to_perturb <=100) system(paste0("echo perturb_interface_4nodes >>",pert))